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Abstract 

Given  an  m  x  n  matrix  A  and  a  positive  integer  k,  we  introduce  a  randomized 
procedure  for  the  approximation  of  A  with  a  matrix  B  of  rank  k.  The  procedure  relies 
on  applying  A"^  to  a  collection  of  I  random  vectors,  where  I  is  an  integer  equal  to 
or  slightly  greater  than  k;  the  scheme  is  efficient  whenever  A  and  A"^  can  be  applied 
rapidly  to  arbitrary  vectors.  The  discrepancy  ||i3  —  A||  is  of  the  same  order  as  the 
{k  +  1)®*  greatest  singular  value  ak+i  of  A,  with  negligible  probability  of  even  moder¬ 
ately  large  deviations.  The  actual  estimates  derived  in  the  paper  are  fairly  complicated, 
but  simplify  when  I  —  k  is  fixed  at  a  small  nonnegative  integer.  For  example,  according 
to  one  of  our  estimates  for  I  —  k  =  20,  the  probability  that  \\B  —  A||  is  greater  than 
10  y^k{k  +  20)  mn  ak+i  is  less  than  10“^^.  The  paper  contains  a  number  of  estimates 
for  II S  — All,  including  several  that  are  stronger  (but  more  detailed)  than  the  preceding 
example.  The  scheme  provides  a  simple,  efficient  means  for  constructing  an  accurate 
approximation  to  the  Singular  Value  Decomposition  of  A,  and  operates  reliably  in¬ 
dependently  of  the  structure  of  the  matrix  A.  The  results  are  illustrated  via  several 
numerical  examples. 


1  Introduction 

In  many  practical  circumstances,  it  is  desirable  to  approximate  a  matrix  A  with  a  sum  of 
rank-1  matrices.  Such  an  approximation  of  A  often  facilitates  understanding  of  the  properties 
of  A.  Moreover,  if  the  approximation  involves  only  a  small  number  of  rank-1  matrices,  then 
the  approximation  also  facilitates  rapid  calculations  involving  A. 

There  are  at  least  two  classical  forms  of  such  matrix  approximations.  One  is  an  ap¬ 
proximation  to  a  Singular  Value  Decomposition  (SVD),  which  is  known  in  the  statistical 
literature  as  a  Principal  Component  Analysis.  The  other  is  an  approximation  obtained  via 
subset  selection;  we  will  refer  to  the  matrix  representation  obtained  via  subset  selection  as 
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Grant  #HR0011-05-l-0002. 
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an  interpolative  decomposition.  These  two  types  of  matrix  approximations  are  defined  as 
follows. 

An  approximation  to  an  SVD  of  a  real  m  x  n  matrix  A  consists  of  nonnegative  real 
numbers  Ui,  <72,  ...,  <Jk-i,  o'k  known  as  singular  values,  orthonormal  real  m  x  1  column 
vectors  u^,  ...,  u^~^,  known  as  left  singular  vectors,  and  orthonormal  real  n  x  1 

column  vectors  . . . ,  known  as  right  singular  vectors,  such  that 


A 


Z 

i=i 


a,- 


(1) 


where  k,  m,  and  n  are  positive  integers,  5  is  a  positive  real  number  specifying  the  precision 
of  the  approximation,  and,  for  any  matrix  B,  ||i?||  denotes  the  spectral  (/^-operator)  norm 
of  B,  that  is,  ||i?||  is  the  greatest  singular  value  of  B.  An  approximation  to  an  SVD  of  A  is 
often  written  in  the  equivalent  form 

\\A-UZV^\\<6,  (2) 


where  U  is  a  real  mxk  matrix  whose  columns  are  orthonormal,  V  is  a  real  nxk  matrix  whose 
columns  are  orthonormal,  and  S  is  a  real  k  x  k  matrix  whose  entries  are  all  nonnegative  and 
whose  entries  off  of  the  main  diagonal  are  zero.  See,  for  example,  [18]  for  a  discussion  of 
SVDs. 

An  interpolative  decomposition  of  a  real  m  x  n  matrix  A  consists  of  a  real  mxk  matrix 
B  whose  columns  constitute  a  subset  of  the  columns  of  A,  and  a  real  k  x  n  matrix  P,  such 
that 

1.  some  subset  of  the  columns  of  P  makes  up  the  k  x  k  identity  matrix, 

2.  no  entry  of  P  has  an  absolute  value  greater  than  2,  and 

3.  A  =  BP. 


See,  for  example,  [17],  [6],  or  Sections  4  and  5  of  [4]  for  a  discussion  of  interpolative  decom¬ 
positions. 

Given  an  algorithm  permitting  the  fast  application  of  a  numerically  low-rank  matrix  A, 
and  an  algorithm  permitting  the  fast  application  of  A"*",  the  algorithm  of  the  present  paper 
provides  a  simple,  efficient  way  for  computing  an  accurate  approximation  to  an  SVD  of  A. 
Moreover,  the  algorithm  provides  a  similar  method  for  computing  an  accurate  approximation 
to  an  interpolative  decomposition  of  A  under  the  same  circumstances. 

Our  scheme  also  provides  an  efficient,  robust  means  for  approximating  the  k  greatest 
singular  values  and  corresponding  singular  vectors  of  any  matrix  A  for  which  a  representation 
enabling  the  fast  application  of  both  A  and  is  available.  The  precision  5  of  the  resulting 
approximation  given  by  formula  (2)  is  at  most  a  reasonably  small  multiple  of  the  {k  + 
l)st  singular  value  of  A.  In  this  regard,  the  algorithm  described  below  should  be 

compared  to  the  classical  Lanczos  method  (for  a  description  of  the  Lanczos  method,  see,  for 
example.  Chapter  9  in  [15]). 

Unlike  the  deterministic  Lanczos  scheme,  the  algorithm  of  the  present  paper  is  a  ran¬ 
domized  one,  and  fails  with  a  rather  negligible  probability.  Examples  of  the  probabilities 
involved  can  be  found  in  (92)  and  (103)  in  Section  4  below. 
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Some  potential  applications  of  the  algorithm  include  hnding  the  eigenmodes  of  certain 
networks,  mining  digital  documents  for  information  via  latent  semantic  analysis,  computing 
electron  densities  within  the  density  functional  theory  of  quantum  chemistry,  simplifying 
the  implementation  of  algorithms  for  fast  matrix  inversion  that  are  based  on  the  compres¬ 
sion  of  blocks  within  matrices,  and  improving  condition  number  estimation  and  subspace 
determination  algorithms  that  are  based  on  inverse  iteration. 

We  should  point  out  that  [16]  and  [4]  motivated  many  aspects  of  the  algorithm  and 
analysis  of  the  present  paper.  Moreover,  a  number  of  recent  publications  address  issues 
similar  to  those  addressed  by  the  present  paper;  we  refer  the  reader  to  [19]  and  [2],  which 
describe  deterministic  methods,  and  to  [1],  [7],  [8],  [9],  [10],  [11],  [12],  [13],  and  the  extensive 
references  contained  therein,  all  of  which  describe  probabilistic  Monte  Carlo  methods. 

We  do  not  analyze  in  detail  the  effects  of  round-off  upon  the  algorithm  of  the  present 
paper.  However,  most  of  the  bounds  that  we  discuss  have  hnite-precision  analogues.  This  is 
conhrmed  by  both  our  preliminary  analysis  and  our  numerical  experiments  (some  of  which 
are  described  in  Section  5  below).  For  simplicity,  we  discuss  only  real  matrices;  the  analysis 
below  extends  easily  to  the  complex  case. 

The  present  paper  has  the  following  structure:  Section  2  collects  together  various  known 
facts  which  later  sections  utilize.  Section  3  provides  the  principal  lemmas  which  Section  4  uses 
to  construct  algorithms.  Section  4  describes  the  algorithm  of  the  present  paper,  providing 
details  about  its  accuracy  and  computational  costs,  and  Section  5  illustrates  the  algorithm 
via  several  numerical  examples. 


2  Preliminaries  from  linear  algebra  and  the  theory  of 
probability 

In  this  section,  we  summarize  various  facts  about  matrices.  Subsection  2.1  discusses  the 
approximation  of  arbitrary  matrices.  Subsection  2.2  discusses  the  singular  values  of  arbitrary 
matrices.  Subsection  2.3  discusses  the  singular  values  of  certain  random  matrices. 

In  the  present  section  and  throughout  the  rest  of  the  paper,  we  employ  the  following 
notation.  In  accordance  with  the  standard  practice,  we  will  denote  the  base  of  the  natural 
logarithm  by  e.  For  any  matrix  A,  we  dehne  the  norm  ||H||  of  A  to  be  the  spectral  (Z^- 
operator)  norm  of  A,  that  is,  ||H||  is  the  greatest  singular  value  of  A.  For  any  positive 
integer  n,  and  real  n  x  1  column  vector  v  G  R”,  we  dehne  the  norm  ||n||  of  v  to  be  the 
root-sum-square  (/^  norm)  of  the  entries  of  v,  that  is. 


T  = 


k=l 


(3) 


where  Vk  is  the  entry  of  v.  (Of  course,  the  norm  of  v  as  viewed  as  a  real  n  x  1  matrix  is 
equal  to  the  norm  of  v  as  viewed  as  a  real  n  x  1  column  vector.) 
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2.1  Approximation  of  general  matrices 

The  following  lemma  states  that,  for  any  m  x  n  matrix  A  whose  rank  is  k,  where  k,  m,  and 
n  are  positive  integers,  there  exist  an  m  x  k  matrix  B  whose  columns  constitute  a  subset  of 
the  columns  of  A,  and  a  k  x  n  matrix  P,  such  that 

1.  some  subset  of  the  columns  of  P  makes  up  the  k  x  k  identity  matrix, 

2.  P  is  not  too  large,  and 

3.  BP  =  A. 

Moreover,  the  lemma  provides  an  analogous  approximation  B  P  to  A  when  the  exact  rank 
of  A  is  not  k,  but  the  {k  +  1)®*  singular  value  of  A  is  nevertheless  small.  The  lemma  is  a 
reformulation  of  Theorem  3.2  in  [17]  and  Theorem  3  in  [6]. 

Lemma  1  Suppose  that  m  and  n  are  positive  integers,  and  A  is  a  real  m  x  n  matrix. 

Then,  for  any  positive  integer  k  with  k  <  m  and  k  <  n,  there  exist  a  real  k  x  n  matrix 
P,  and  a  real  m  x  k  matrix  B  whose  columns  constitute  a  subset  of  the  columns  of  A,  such 
that 

1.  some  subset  of  the  columns  of  P  makes  up  the  k  x  k  identity  matrix, 

2.  no  entry  of  P  has  an  absolute  value  greater  than  1, 

3.  ||P||  <  ^/k{n-  k)  +  1, 

4-  the  least  (that  is,  the  k^^  greatest)  singular  value  of  P  is  at  least  1, 

5.  B  P  =  A  when  k  =  m  or  k  =  n,  and 

6.  ||P  P  —  A||  <  ^Jk{n  —  k)  +  1  (jfc+i  when  k  <  m  and  k  <  n,  where  is  the  {k  +  1)^* 
greatest  singular  value  of  A. 

Remark  2  Properties  1,  2,  3,  and  4  in  Lemma  1  ensure  that  the  interpolative  decomposition 
P  P  of  A  is  numerically  stable.  Also,  Property  3  follows  directly  from  Properties  1  and  2, 
and  Property  4  follows  directly  from  Property  1. 

Observation  3  There  exists  an  algorithm  which  computes  B  and  P  in  Lemma  1  from  A, 
provided  that  we  require  only  that 

1.  some  subset  of  the  columns  of  P  makes  up  the  k  x  k  identity  matrix, 

2.  no  entry  of  P  has  an  absolute  value  greater  than  2, 

3.  ||P||  <  A/4F(n'^^X)”TT, 

4.  the  least  (that  is,  the  greatest)  singular  value  of  P  is  at  least  1, 

5.  B  P  =  A  when  k  =  m  or  k  =  n,  and 
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6.  \\B  P  —  A||  <  ^jAk{n  —  k)  +  1  cifc+i  when  k  <  m  and  k  <  n,  where  ak+i  is  the  {k  +  1)®* 
greatest  singular  value  of  A. 

For  any  positive  real  number  £,  the  algorithm  can  identify  the  least  k  such  that  ||-B  -P— ^||  ~  £• 
Furthermore,  there  exists  a  real  number  C  such  that  the  algorithm  computes  both  B  and 
P  using  at  most  Ckmn  log(min(m,  n))  floating-point  operations  and  Ckmn  floating-point 
words  of  memory.  The  algorithm  is  based  upon  the  Cramer  rule  and  the  ability  to  obtain 
the  minimal-norm  (or  at  least  roughly  minimal-norm)  solutions  to  linear  algebraic  systems 
of  equations  (see  [17],  [6],  and  [16]). 


2.2  Singular  values  of  general  matrices 

The  following  technical  lemma  will  be  needed  in  Section  3. 

Lemma  4  Suppose  that  m  and  n  are  positive  integers  with  m  >  n.  Suppose  further  that  A 
is  a  real  m  x  n  matrix  such  that  AA  A  is  invertible. 

Then, 

t|(.4T.4)-‘AT||  =  -,  (4) 

where  an  is  the  least  (that  is,  the  greatest)  singular  value  of  A. 

Proof.  We  form  an  SVD  of  A, 

A  =  U^V^,  (5) 

where  7/  is  a  real  unitary  mxm  matrix,  S  is  a  real  mxn  matrix  whose  entries  are  nonnegative 
everywhere  and  zero  off  of  the  main  diagonal,  and  P  is  a  real  unitary  nxn  matrix.  It  follows 
from  (5)  that 

{A^  A)-^  A^  =  v  (6) 

Combining  (6)  and  the  fact  that  U  and  V  are  unitary  yields 

||(7l^7l)-M^||  =  ||(S^S)-^S^||  .  (7) 

Combining  (7)  and  the  fact  that  S  is  zero  off  of  the  main  diagonal  yields  (4).  □ 

The  following  lemma  provides  what  is  known  as  the  Courant-Fischer  maximin  character¬ 
ization  of  singular  values;  Theorem  8.1.2  in  [15]  provides  an  equivalent  formulation  of  (8). 


Lemma  5  Suppose  that  m  and  n  are  positive  integers,  and  A  is  a  real  mxn  matrix. 
Then,  the  k^^  greatest  singular  value  ak  of  A  is  given  by  the  formula 


\\Av\ 

ak  =  max  mm  — ^ 

SCR":  dim 5=^  -uSS:  ||ti||7^0  U 


for  all  k  =  1,  2,  ...,  min(m,n)  —  1,  min(m, n),  where  the  maximum  is  taken  over  all  k- 
dimensional  subspaces  of  E”,  and  the  minimum  is  taken  over  all  vectors  in  S  that  have 
nonzero  norms. 
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The  following  lemma  states  that  the  singular  values  of  the  product  G  A  of  matrices  G 
and  A  are  at  most  ||G||  times  greater  than  the  corresponding  singular  values  of  A. 


Lemma  6  Suppose  that  I,  m,  and  n  are  positive  integers,  A  is  a  real  m  x  n  matrix,  and  G 
is  a  real  I  x  m  matrix. 

Then,  the  greatest  singular  value  pk  of  the  product  G  A  is  at  most  a  factor  of  ||G|| 
times  the  greatest  singular  value  of  A,  that  is, 

Pk  <  ||G||  O'k  (9) 

for  all  k  =  1,  2,  . . . ,  min(/,  m,  n)  —  1,  min(/,  m,  n) . 


Proof.  For  any  vector  v  G  R"'  with  ||r;||  ^  0, 


(10) 


Combining  (8)  and  (10)  yields  (9). 


□ 


The  following  lemma  states  that  the  greatest  singular  value  of  a  matrix  A  is  at  least  as 
large  as  the  greatest  singular  value  of  any  rectangular  block  of  entries  in  A]  the  lemma  is  a 
straightforward  consequence  of  the  minimax  properties  of  singular  values  (see,  for  example, 
Section  47  of  Chapter  2  in  [20]). 


Lemma  7  Suppose  that  k,  I,  m,  and  n  are  positive  integers  with  k  <m  and  I  <  n.  Suppose 
further  that  A  is  a  real  m  x  n  matrix,  and  B  is  a  k  x  I  rectangular  block  of  entries  in  A. 
Then,  the  greatest  singular  value  of  B  is  at  most  the  greatest  singular  value  of  A. 


The  following  lemma  states  that  if  the  norm  of  the  difference  of  two  matrices  is  small,  then 
their  corresponding  singular  values  are  close;  Corollary  8.6.2  in  [15]  provides  an  equivalent 
formulation  of  (11). 

Lemma  8  Suppose  that  m  and  n  are  positive  integers,  and  A  and  R  are  real  mxn  matrices. 

Then,  the  k^^  greatest  singular  value  Tk  of  the  sum  A  +  R  and  the  k^^  greatest  singular 
value  at  of  A  differ  by  at  most  ||i?||,  that  is, 

\rk  -  cTkl  <  \\R\\  (11) 

for  all  k  =  1,  2,  . . . ,  min(m,  n)  —  1,  min(m,  n) . 


2.3  Singular  values  of  random  matrices 

The  following  lemma  provides  a  highly  probable  upper  bound  on  the  greatest  singular  value 
of  a  square  matrix  whose  entries  are  independent,  identically  distributed  (i.i.d.)  Gaussian 
random  variables  of  zero  mean  and  unit  variance;  Formula  8.8  of  [14]  provides  an  equivalent 
formulation  of  the  lemma. 
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Lemma  9  Suppose  that  n  is  a  positive  integer,  G  is  a  real  n  x  n  matrix  whose  entries  are 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  varianee,  and  j  is  a  positive  real 
number,  such  that  7  >  1  and 

272 

£7^-1 

is  nonnegative. 

Then,  the  greatest  singular  value  of  G  is  at  most  v^^7  with  probability  not  less  than  the 
amount  in  (12). 

Combining  Lemmas  7  and  9  yields  the  following  lemma,  providing  a  highly  probable  upper 
bound  on  the  greatest  singular  value  of  a  rectangular  matrix  whose  entries  are  i.i.d.  Gaussian 
random  variables  of  zero  mean  and  unit  variance. 

Lemma  10  Suppose  that  I,  m,  and  n  are  positive  integers  with  n  >  I  and  n  >  m.  Suppose 
further  that  G  is  a  real  I  x  m  matrix  whose  entries  are  i.i.d.  Gaussian  random  variables  of 
zero  mean  and  unit  variance,  and  ■y  is  a  positive  real  number,  such  that  7  >  1  and  (12)  is 
nonnegative. 

Then,  the  greatest  singular  value  of  G  is  at  most  v^^7  with  probability  not  less  than  the 
amount  in  (12). 

The  following  lemma  provides  a  highly  probable  lower  bound  on  the  least  singular  value 
of  a  rectangular  matrix  whose  entries  are  i.i.d.  Gaussian  random  variables  of  zero  mean  and 
unit  variance;  Formula  2.5  in  [5]  and  the  proof  of  Lemma  4.1  in  [5]  together  provide  an 
equivalent  formulation  of  Lemma  11. 

Lemma  11  Suppose  that  k  and  I  are  positive  integers  with  k  <  1.  Suppose  further  that  G 
is  a  real  I  x  k  matrix  whose  entries  are  i.i.d.  Gaussian  random  variables  of  zero  mean  and 
unit  variance,  and  (3  is  a  positive  real  number,  such  that 

(l-k  +  i))]) 

is  nonnegative. 

Then,  the  least  (that  is,  the  k*^  greatest)  singular  value  of  G  is  at  least  /3)  with 

probability  not  less  than  the  amount  in  (13). 

3  Mathematical  apparatus 

In  this  section,  we  describe  the  principal  tools  used  in  Section  4. 

The  following  lemma  states  that  the  product  B  P  oi  matrices  B  and  P  is  a  good  approx¬ 
imation  to  a  matrix  A,  provided  that  there  exists  a  matrix  G  such  that 

1.  the  columns  of  B  constitute  a  subset  of  the  columns  of  A, 

2.  ||P||  is  not  too  large. 


1  - 


4  (72  —  1) 
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3.  G5P  is  a  good  approximation  to  G  A,  and 

4.  there  exists  a  matrix  F  such  that  ||P||  is  not  too  large,  and  PGA  is  a  good  approxi¬ 
mation  to  A. 

Lemma  12  Suppose  that  k,  I,  m,  and  n  are  positive  integers  with  k  <  n.  Suppose  further 
that  A  is  a  real  m  x  n  matrix,  B  is  a  real  mxk  matrix  whose  columns  constitute  a  subset  of 
the  columns  of  A,  P  is  a  real  k  x  n  matrix,  F  is  a  real  m  x  I  matrix,  and  G  is  a  real  I  x  m 
matrix. 

Then, 

\\BP-A\\  <  \\FGA-A\\{\\P\\  +  1)  +  \\F\\\\GBP-GA\\.  (14) 

Proof.  We  observe  that 

IIP  P- All  <  \\B  P  -  F G  B  P\\  +  \\F  G  B  P  -  F  G  A\\  +  \\F  G  A  -  A\\,  (15) 

\\BP-FGBP\\<\\B-FGB\\\\P\\,  (16) 

and 

\\FGBP-FGA\\  <  ||P||  \\GBP-GA\\.  (17) 

Since  the  columns  of  P  constitute  a  subset  of  the  columns  of  A,  it  follows  that  the  columns 
of  B  —  F  G  B  constitute  a  subset  of  the  columns  of  A  —  P  G  A,  and  therefore, 

\\B-FGB\\<\\A-FGA\\.  (18) 

Combining  (15),  (16),  (17),  and  (18)  yields  (14).  □ 


Remark  13  Since  the  columns  of  P  constitute  a  subset  of  the  columns  of  A  in  Lemma  12, 
it  follows  that  the  columns  of  GP  constitute  a  subset  of  the  columns  of  GA.  Conversely, 
whenever  a  matrix  S  is  formed  by  gathering  distinct  columns  of  GA  together  into  S,  then 
clearly  S  =  G  B  for  some  matrix  P  whose  columns  constitute  a  subset  of  the  columns  of  A. 

The  following  lemma  states  that,  for  any  matrix  A,  and  matrix  G  whose  entries  are 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  with  very  high  probability 
there  exists  a  matrix  P  with  a  reasonably  small  norm,  such  that  PGA  is  a  good  approxi¬ 
mation  to  A. 


Lemma  14  Suppose  that  k,  I,  m,  and  n  are  positive  integers  with  k  <  I,  such  that  I  <  m 
and  I  <  n.  Suppose  further  that  A  is  a  real  m  x  n  matrix,  G  is  a  real  I  x  m  matrix  whose 
entries  are  i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  and  f3  and  7 
are  positive  real  numbers,  such  that  7  >  1  and 


^271  {I  —  k  +  1)  \{l  —  k  +  1)  (3 


4  (7^  —  1)  \G 


is  nonnegative. 


Then,  there  exists  a  real  m  x  I  matrix  F  such  that 


\\FGA-A\\  <  V2/m/32 72  +  1  afc+i  (20) 

and 

\\F\\<Vifi  (21) 

with  probability  not  less  than  the  amount  in  (19),  where  cifc+i  is  the  {k  +  1)®*  greatest  singular 
value  of  A. 

Proof.  We  prove  the  existence  of  a  matrix  F  satisfying  (20)  and  (21)  by  constructing  one. 
We  start  by  forming  an  SVD  of  A, 


A  =  UTV^,  (22) 

where  is  a  real  unitary  mxm  matrix,  S  is  a  real  mxn  matrix  whose  entries  are  nonnegative 
everywhere  and  zero  off  of  the  main  diagonal,  and  P  is  a  real  unitary  n  x  n  matrix,  such 
that 

(23) 

for  alH  =  1,  2,  . . . ,  min(m,n)  —  1,  min(m,n),  where  Sj  j  is  the  entry  in  row  i  and  column  i 
of  S,  and  a*  is  the  greatest  singular  value  of  A. 

Next,  we  dehne  auxiliary  matrices  H,  R,  and  P.  We  dehne  H  to  be  the  leftmost  /  x  k 
block  of  the  I  x  m  matrix  G  U,  and  R  to  be  the  rightmost  I  x  [m  —  k)  block  of  G  U,  so  that 

GU  =  {  H  \  R)  .  (24) 

Combining  the  facts  that  U  is  real  and  unitary,  and  that  the  entries  of  G  are  i.i.d.  Gaussian 
random  variables  of  zero  mean  and  unit  variance,  we  see  that  the  entries  of  H  are  also 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  as  are  the  entries  of  R.  We 
dehne  to  be  the  real  k  x  I  matrix  given  by  the  formula 

=  {H^  H)-^  .  (25) 

We  dehne  P  to  be  the  m  x  I  matrix  whose  uppermost  k  x  I  block  is  H^~^\  and  whose  entries 
in  the  lowermost  {m  —  k)  x  I  block  are  all  zeros,  so  that 

P  =  (^)  ■  (26) 

Finally,  we  dehne  F  to  be  the  m  x  I  matrix  given  by 

F  =  UF^U  (27) 

Combining  (25),  (4),  the  fact  that  the  entries  of  H  are  i.i.d.  Gaussian  random  variables 
of  zero  mean  and  unit  variance,  and  Lemma  11  yields 

fd  (28) 
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with  probability  not  less  than 


1 


1 

\/2tt  {I  —  k  +  1) 


k  +  l)/]) 


(29) 


Combining  (27),  (28),  and  the  fact  that  U  is  unitary  yields  (21). 

We  now  show  that  F  dehned  in  (27)  satishes  (20). 

We  dehne  S  to  be  the  leftmost  uppermost  k  x  k  block  of  S,  and  T  to  be  the  rightmost 
lowermost  {m  —  k)  x  [n  —  k)  block  of  S,  so  that 


S  = 


Combining  (22),  (24),  and  (27)  yields 
FGA-A=U 

Combining  (25)  and  (30)  yields 

//(-!) 


'  .s 

0 

0 

T 

-1)  ' 

L 

0 


0 


H  R 

)-lj  SC 

H^-^^RT 

VO 

-T 

Furthermore, 


0 

R  T  ' 

0 

-T 

<  RT  ^ +  \\Tf 


Moreover, 

Combining  (30)  and  (23)  yields 

l|r||  <  fft+i. 

Combining  (31),  (32),  (33),  (34),  (35),  and  the  fact  that  U  and  V  are  unitary  yields 


(30) 

(31) 

(32) 

(33) 

(34) 

(35) 


||fG/l-/l||  <  llfif  +  lot+i. 


(36) 


Combining  Lemma  10  and  the  fact  that  the  entries  of  R  are  i.i.d.  Gaussian  random 
variables  of  zero  mean  and  unit  variance  shows  that 


\\R\\  <  7 


with  probability  not  less  than 

1 _ ^ _ (2A-) 

4  (72  —  1)  ^  J 


Combining  (36),  (28),  and  (37)  yields  (20). 


(37) 


(38) 

□ 


The  following  lemma  is  very  similar  to  Lemma  14.  Lemma  15  is  tighter  than  Lemma  14 
when  the  singular  values  of  the  matrix  A  decay  sufficiently  fast,  and  the  numbers  j  and  /  in 
the  lemma  are  both  much  less  than  m. 
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Lemma  15  Suppose  that  j,  k,  I,  m,  and  n  are  positive  integers  with  k  <  I,  such  that 
k  +  j  <  m  and  k  +  j  <  n,  as  well  as  I  <  m  and  I  <  n.  Suppose  further  that  A  is  a  real  mxn 
matrix,  G  is  a  real  I  x  m  matrix  whose  entries  are  i.i.d.  Gaussian  random  variables  of  zero 
mean  and  unit  variance,  and  f3  and  7  are  positive  real  numbers,  such  that  7  >  1  and 


$  =  1  - 


y^27r  {I  —  k  +  1)  \{l  —  k  +  1)  (3 

1 


4  (7^  —  1)  a/tt  max(m  —  k  —  1)  7^ 


27^ 


2  \  max(m— fc— 


2-1 


27^ 

4  (7^  —  1)  A/7r^nax(70^ 


2  \  max(j,  Z) 


(39) 


is  nonnegative. 

Then,  there  exists  a  real  m  x  I  matrix  F  such  that 


F  G  A  — A\\  <  a/2/  max(j,  /)  72  -|-  i  +  \/2l  max(m  —  k  —  j,  1)  7^  +  1  (40) 


and 

\\F\\<Vifi  (41) 

with  probability  not  less  than  the  amount  in  (39),  where  cifc+i  is  the  {k  +  1)®*  greatest  singular 
value  of  A,  and  o-fc+j+i  *■5  the  {k  +  j  +  1)^*  greatest  singular  value  of  A. 


Proof.  We  prove  the  existence  of  a  matrix  F  satisfying  (40)  and  (41)  by  constructing  one. 
We  start  by  forming  an  SVD  of  A, 


A  =  UTV^, 


(42) 


where  U  is  a  real  unitary  mxm  matrix,  S  is  a  real  mxn  matrix  whose  entries  are  nonnegative 
everywhere  and  zero  off  of  the  main  diagonal,  and  P  is  a  real  unitary  n  x  n  matrix,  such 
that 

=  ai  (43) 

for  alH  =  1,  2,  . . . ,  min(m,n)  —  1,  min(m,n),  where  Sj  j  is  the  entry  in  row  i  and  column  i 
of  S,  and  cij  is  the  greatest  singular  value  of  A. 

Next,  we  dehne  auxiliary  matrices  H,  R,  F,  and  P.  We  dehne  H  to  be  the  leftmost  /  x  k 
block  of  the  /  x  m  matrix  GU ,  R  io  be  the  I  x  j  block  of  GU  whose  hrst  column  is  the 
{k  +  1)®*  column  of  G  U,  and  F  to  be  the  rightmost  lx  {m  —  j  —  k)  block  of  G  U,  so  that 

GU  =  {  H  \  R\T  )  .  (44) 

Combining  the  facts  that  U  is  real  and  unitary,  and  that  the  entries  of  G  are  i.i.d.  Gaussian 
random  variables  of  zero  mean  and  unit  variance,  we  see  that  the  entries  of  H  are  also 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  as  are  the  entries  of  R,  and 
as  are  the  entries  of  F.  We  dehne  to  be  the  real  k  x  I  matrix  given  by  the  formula 

=  {H^  H)-^  .  (45) 
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We  define  P  to  be  the  real  m  x  /  matrix  whose  uppermost  k  x  /  block  is  H^~^\  whose  entries 
in  the  j  x  /  block  whose  hrst  row  is  the  {k  +  1)®*  row  of  P  are  all  zeros,  and  whose  entries  in 
the  lowermost  {m  —  k  —  j)  x  I  block  are  all  zeros,  so  that 


P  = 


0 


(46) 


Finally,  we  dehne  F  to  be  the  m  x  I  matrix  given  by 


F  =  UP 


U 


0 

0 


(47) 


Combining  (45),  (4),  the  fact  that  the  entries  of  H  are  i.i.d.  Gaussian  random  variables 
of  zero  mean  and  unit  variance,  and  Lemma  11  yields 


<Vlp 


with  probability  not  less  than 

1 


\/2tt  {I  —  k  +  1)  \{l  —  k  +  1)  (3 


(48) 


(49) 


Combining  (47),  (48),  and  the  fact  that  U  is  unitary  yields  (41). 

We  now  show  that  F  dehned  in  (47)  satisfies  (40). 

We  dehne  S  to  be  the  leftmost  uppermost  k  x  k  block  of  S,  T  to  be  the  j  x  j  block  of  S 
whose  leftmost  uppermost  entry  is  the  entry  in  the  {k  +  1)®*  row  and  {k  +  1)®*  column  of  S, 
and  0  to  be  the  rightmost  lowermost  {m  —  k  —  j)  x  {n  —  k  —  j)  block  of  S,  so  that 


S  = 


'  5 

0 

0  ' 

0 

T 

0 

0 

0 

0 

(50) 


Combining  (42),  (44),  and  (47)  yields 
FGA-A=U 

Combining  (45)  and  (50)  yields 

'  (  ^  I  ^  I  r  ) 


0 

6 


H  \R\T  )  -1  S  1/^ 


(51) 


0 

0 


s  = 


0 

H^-^^RT 

i7(-i)  r  0  ' 

0 

-T 

0 

0 

0 

-0 

(52) 


Furthermore, 


0 

H^-^^RT 

//(-i)  r  0  ' 

0 

-T 

0 

0 

0 

-0 

<  F  0  V  ||T||2  +  ||0||2.  (53) 


12 


Moreover, 

||/?i|  ||T||  (54) 

and 

r  0||  <  ||r||  ||0||.  (55) 

Combining  (50)  and  (43)  yields 

ni  <  ffs+i  (56) 

and 

lien  <  <^<=+,+1-  (5T) 

Combining  (51),  (52),  (53),  (54),  (55),  (56),  (57),  and  the  fact  that  U  and  V  are  unitary 
yields 

\\FGA-Af<  Pir  +  l)  (afc+i)2+ ||rf  +  l)  {ak+j+if.  (58) 

Combining  Lemma  10  and  the  fact  that  the  entries  of  R  are  i.i.d.  Gaussian  random 
variables  of  zero  mean  and  unit  variance,  as  are  the  entries  of  L,  yields 

\\R\\  <  ^2  max(j,/)  7  (59) 

and 

||r||  <  a/2  max(m  —  k  —  j,  1)  7,  (60) 

with  probability  not  less  than 


27^ 

4  (7^  —  1)  a/tt  max(m  —  k  —  j,  1)  7^ 


2  \  max(m— /c— Z) 


27^ 

4  (72  —  1)  i/tt  max(j,  /)  7^ 
Combining  (58),  (48),  (59),  and  (60)  shows  that 


2  \  max(7Z) 


•  (61) 


\\FG  A-A\\'^  <  (2/  max(j,/)/3^7^  +  1)  {ak+if  +  {2l  max(m  - /c  -  j, /) 7^  +  l)  {ak+j+if 

(62) 

with  probability  not  less  than  the  amount  in  (39).  Combining  (62)  and  the  fact  that 


\/x  +  y  <  y/x  +  y/y 

for  any  nonnegative  real  numbers  x  and  y  yields  (40). 


(63) 

□ 


Given  an  m  x  n  matrix  A,  and  a  matrix  G  whose  entries  are  i.i.d.  Gaussian  random 
variables  of  zero  mean  and  unit  variance,  the  following  lemma  provides  an  upper  bound  on 
the  singular  values  of  the  product  GA  in  terms  of  the  singular  values  of  A]  the  lemma  is 
most  useful  when  the  singular  values  of  A  decay  sufficiently  fast,  and  the  numbers  j  and  I 
in  the  lemma  are  both  much  less  than  m. 
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Lemma  16  Suppose  that  j,  k,  I,  m,  and  n  are  positive  integers  with  k  <  I,  such  that 
k  +  j  <  m  and  k  +  j  <  n.  Suppose  further  that  A  is  a  real  m  x  n  matrix,  G  is  a  real  I  x  m 
matrix  whose  entries  are  i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance, 
and  'J  is  a  positive  real  number,  such  that  7  >  1  and 

-1  /  o  9  \  max(m— fc— 7,  0 

_ 1  \ 

4  (72  —  1)  a/ tt  max(m  —  k  —  j,  1)  7^  ^  / 

_ 1 _  f  27^ 

4  (72  —  1)  a/tt  max(A;  +  j,  /)  7^ 


max(/c+j',  /) 


(64) 


is  nonnegative. 

Then,  the  {k  +  1)®*  greatest  singular  value  pk+i  of  G  A  is  at  most  a  certain  linear  combi¬ 
nation  of  the  {k  +  1)®*  greatest  singular  value  ak+i  of  A  and  the  {k-\-j-\- 1)®*  greatest  singular 
value  Ck+j+i  of  A,  namely. 


Pk+i  <  a/2  max(/c  +  j,  /)  7  au+i  +  a/2  max(m  -  k-  j,l)  'y  (Jk+j+i,  (65) 

with  probability  not  less  than  the  amount  in  (64)- 

Proof.  We  start  by  forming  an  SVD  of  A, 

A  =  Ui:V^,  (66) 

where  U  is  a  real  unitary  mxm  matrix,  S  is  a  real  mxn  matrix  whose  entries  are  nonnegative 
everywhere  and  zero  off  of  the  main  diagonal,  and  P  is  a  real  unitary  n  x  n  matrix,  such 
that 

=  ai  (67) 

for  alH  =  1,  2,  . . . ,  min(m,n)  —  1,  min(m,n),  where  is  the  entry  in  row  i  and  column  i 
of  S,  and  a*  is  the  greatest  singular  value  of  A. 

Combining  (66)  and  the  fact  that  V  is  unitary  yields  that  GA  has  the  same  singular 
values  as  GUT,. 

Next,  we  dehne  auxiliary  matrices  H  and  R.  We  dehne  H  to  be  the  leftmost  lx{k  +  j) 
block  of  the  /  x  m  matrix  G  U,  and  R  to  be  the  rightmost  I  x  [m  —  k  —  j)  block  of  G  U,  so 
that 

GU={H\R).  (68) 

Combining  the  facts  that  U  is  real  and  unitary,  and  that  the  entries  of  G  are  i.i.d.  Gaussian 
random  variables  of  zero  mean  and  unit  variance,  we  see  that  the  entries  of  H  are  also 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  as  are  the  entries  of  R. 

Gombining  (68)  and  the  fact  that  G  A  has  the  same  singular  values  as  GUT  yields  that 
G  A  has  the  same  singular  values  as  (  77  |  0  )  S+(0|7?)  S. 

It  follows  from  (67)  that 


(  0  I  7?  )  S||  <  ||77||  ak+j+i. 


(69) 
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Combining  (11)  and  (69)  yields 


Pk+l  <  'Tk+l  +  ||-R||  (^k+j+l)  (70) 

where  p^+i  is  the  {k  +  1)®*  greatest  singular  value  of(if|0)S  +  (0|-R)S,  and  r^+i  is 
the  (A;  + 1)®*  greatest  singular  value  of  (  77  |  0  )  S;  p^+i  is  also  the  {k  +  1)*^*  greatest  singular 
value  of  G  since  G  A  has  the  same  singular  values  as(i7|0)  S+(0|i?)  S. 
Furthermore, 

IK  |0  )||  <  ||ff||.  (71) 

Combining  (9),  (71),  and  (67)  yields 

T-fc+i  <  ||77||  o-fc+i.  (72) 

Combining  Lemma  10  and  the  fact  that  the  entries  of  H  are  i.i.d.  Gaussian  random 
variables  of  zero  mean  and  unit  variance,  as  are  the  entries  of  7?,  shows  that 

||77||  <  a/2  max(/i;  +  j,/)  7  (73) 

and 

||7?||  <  a/2  max(m  —  k  —  j,  1)  7,  (74) 

with  probability  not  less  than  the  amount  in  (64). 

Combining  (70),  (72),  (73),  and  (74)  yields  (65).  □ 

The  following  lemma  provides  an  efficient  means  of  computing  an  SVD  of  a  real  m  x  n 
matrix  A  from  a  real  m  x  k  matrix  B  and  a  real  k  x  n  matrix  P  such  that  A  =  B  P  and  k 
is  much  less  than  both  m  and  n.  If,  in  addition,  ||77||  <  ||A||  and  ||P||  is  not  too  large,  then 
the  scheme  described  by  the  lemma  is  numerically  stable. 


Lemma  17  Suppose  that  k,  m,  and  n  are  positive  integers  with  k  <n.  Suppose  further  that 
A  is  a  real  m  x  n  matrix,  B  is  a  real  m  x  k  matrix,  and  P  is  a  real  k  x  n  matrix,  such  that 


A  =  BP. 


(75) 


Suppose  in  addition  that  L  is  a  real  kxk  matrix,  and  Q  is  a  real  nxk  matrix  whose  columns 
are  orthonormal,  such  that 

P  =  LQ^.  (76) 


Suppose  finally  that  G  is  a  real  m  x  k  matrix,  U  is  a  real  m  x  k  matrix  whose  columns  are 
orthonormal,  T,  is  a  real  kxk  matrix,  and  W  is  a  real  kxk  matrix  whose  columns  are 
orthonormal,  such  that 


and 


=  BL 

(77) 

UEW^. 

(78) 

(79) 
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Then, 


where  V  is  the  real  n  x  k  matrix  given  by  the  formula 


V  =  QW.  (80) 

Moreover,  the  eolumns  ofV  are  orthonormal  (as  are  the  columns  ofU),  and 

l|i|l  =  lli’ll-  (81) 

Proof.  Combining  (75),  (76),  (77),  (78),  and  (80)  yields  (79).  Combining  (80)  and  the 
facts  that  W  is  unitary  and  the  columns  of  Q  are  orthonormal  yields  that  the  columns  of 
V  are  orthonormal.  Combining  (76)  and  the  fact  that  the  columns  of  Q  are  orthonormal 

yields  (81).  □ 


Remark  18  The  matrices  L  and  Q  in  (76)  can  be  computed  from  P  as  follows.  Using  the 
algorithms  described,  for  example,  in  Chapter  5  of  [15],  we  construct  an  upper  triangular 
real  k  x  k  matrix  R,  and  a  real  n  x  k  matrix  Q  whose  columns  are  orthonormal,  such  that 

=  QR.  (82) 

We  thus  obtain  Q.  We  then  dehne  L  to  be  the  transpose  of  R,  that  is, 

L  =  R^.  (83) 

4  Description  of  the  algorithm 

In  this  section,  we  describe  the  algorithm  of  the  present  paper.  In  Subsection  4.1,  we  discuss 
approximations  to  interpolative  decompositions.  In  Subsection  4.2,  we  discuss  approxima¬ 
tions  to  SVDs.  In  Subsection  4.3,  we  tabulate  the  computational  costs  of  various  parts  of 
the  algorithm.  In  Subsection  4.4,  we  describe  Table  1. 

4.1  Interpolative  decomposition 

Suppose  that  k,  m,  and  n  are  positive  integers  with  k  <m  and  k  <  n,  and  4  is  a  real  mxn 
matrix.  In  this  subsection,  we  will  collect  together  k  appropriately  chosen  columns  of  A  into 
a  real  m  x  k  matrix  B,  and  construct  a  real  k  x  n  matrix  P,  such  that 

||P||  <  ^Ak  {n-k)  +  l  (84) 

and 

\\B  P  -  A\\  <  (Jk+i,  (85) 

where  is  the  {k  +  1)®*  greatest  singular  value  of  A.  To  do  so,  we  select  an  integer  I  with 
I  >  k,  such  that  I  <  m  and  /  <  n  (/  =  /c  -|-  20  is  often  a  suitable  choice),  and  make  the 
following  three  steps: 
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1.  Using  a  random  number  generator,  form  a  real  I  x  m  matrix  G  whose  entries  are 
i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  and  compute  the 
I  X  n  product  matrix 

R  =  GA.  (86) 

2.  Using  the  algorithm  of  [6],  form  a  real  I  x  k  matrix  S  whose  columns  constitute  a  subset 
of  the  columns  of  R,  and  a  real  k  x  n  matrix  P  satisfying  (84),  such  that 


||SP-fl||  <  s/ik(n-k)  +  lpM, 


(87) 


where  pk+i  is  the  {k  +  1)*^*  greatest  singular  value  of  R.  (See  Observation  3  for  a  brief 
discussion  of  the  properties  of  the  algorithm  of  [6].) 

3.  Using  the  fact  that  the  columns  of  S  constitute  a  subset  of  the  columns  of  R,  for  any 
j  =  1,  2,  . . . ,  k  —  1,  k,  let  ij  denote  an  integer  such  that  the  column  of  S  is  the 
column  of  R.  Form  the  real  m  x  k  matrix  B  whose  column  is  the  column  of  A 
for  all  j  =  1,  2,  . . . ,  A;  —  1,  k. 

The  matrices  B  and  P  obtained  via  the  preceding  three  steps  satisfy  (84)  and  (85);  see  the 
following  remark. 

Remark  19  In  this  remark,  we  demonstrate  that  the  matrices  B  and  P  satisfy  (85).  Indeed, 
combining  (86)  and  Remark  13  yields 


S  =  GB. 

Combining  (87),  (86),  and  (88)  yields 

IIGBP-GAII  <  s/4k(n-k)  +  lpM, 


(88) 


(89) 


where  pk+i  is  the  {k  +  1)®*  greatest  singular  value  of  R.  Suppose  that  f3  and  7  are  positive 
real  numbers  such  that  7  >  1  and 


X  =  1 


kj2'K  [I  —  k  +  1)  \{l  —  k  +  1)  (3 


l-k+l 


27^ 

2  (72  —  1) 


(90) 


is  nonnegative.  Then,  combining  (14),  (20),  (21),  (84),  (89),  (9),  (86),  and  Lemma  10  yields 

\\BP- A\\ 

<  (^\/2lm(3‘^  7^  +  1  (y\/^k  {n  —  k)  +  1  +  1  j  +  /?  7  V2lm  \/4:k{n  —  k)  +  1  j  (91) 

with  probability  not  less  than  y  defined  in  (90),  where  ak+i  is  the  {k  +  1)®*  greatest  singular 
value  of  A.  The  bound  (91)  is  a  precise  version  of  (85).  For  example,  choosing  (3  =  3/4, 
7^  =  5,  and  /  =  /c  +  20,  and  combining  (91)  and  (90),  we  obtain 


\\B  P  —  A\\  <  10  \/k  {k  +  20)  mn  ak+i 


(92) 


with  probability  not  less  than  1  —  10  Table  1  contains  similar  results  obtained  by  taking 
other  values  for  I  —  k,  P,  and  7. 
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Remark  20  When  the  singular  values  of  A  decay  sufficiently  fast,  and  I  is  much  less  than 
m,  the  factors  ^/2lm(3'^  7^  +  1  and  y/2lm  in  (91)  are  much  larger  than  necessary.  Indeed, 
suppose  that  j  is  a  positive  integer  with  k  +  j  <  m  and  k  +  j  <  n,  and  f3  and  7  are  positive 
real  numbers,  such  that  7  >  1  and  <l>  +  'h  >  1,  where  <h  is  dehned  in  (39),  and  is  dehned 
in  (64).  Then,  combining  (14),  (40),  (41),  (84),  (89),  (65),  and  (86)  yields 

\\BP  -  A\\<^ak+i+r]ak+j+i  (93) 

with  probability  not  less  than  $  +  T  —  1,  where  <h  is  dehned  in  (39),  T  is  dehned  in  (64), 
(Tfc+i  is  the  {k  +  1)®*  greatest  singular  value  of  A,  and  a^+j+i  is  fhe  {k  +  j  +  1)®*  greatest 
singular  value  of  A,  and  where 

^  =  a/2/  max(j,  1)  72  -|-  1  Ak  {n  —  k)  +  1  +  ij 

+  /5  7  a/2/  max(A;  +  j,  /)  a/4A;  (n  —  /c)  +  1  (94) 

and 

7  =  a/2/  max(m  —  k  —  j,  1)  /5^7^  +  1  ^a/4A;  (n  —  /c)  +  1  +  1  j 

+  /97v^  max(m  —  k  -  j,  /)  v^4A;  {n-k)  +  1.  (95) 

When  j,  k,  and  /  are  all  much  less  than  m,  clearly  ^  is  much  less  than  rj. 

Remark  21  If  we  choose  I  =  k  in  the  algorithm  of  the  present  subsection  (instead  of 
choosing  I  >  k),  then  we  must  replace  (87)  with  the  formula 

||^P-R||=0.  (96) 

All  other  aspects  of  the  algorithm  stay  the  same  in  the  case  that  I  =  k.  In  particular,  (91) 
and  (93)  hold  in  the  case  that  /  =  k,  too. 

4.2  Singular  Value  Decomposition 

Suppose  that  k,  m,  and  n  are  positive  integers  with  k  <m  and  k  <  n,  and  A  is  a  real  mxn 
matrix.  In  this  subsection,  we  will  compute  an  approximation  to  an  SVD  of  A  such  that 

\\UJ:V^-A\\<ak+i,  (97) 

where  /7  is  a  real  m  x  k  matrix  whose  columns  are  orthonormal,  is  a  real  n  x  k  matrix 
whose  columns  are  orthonormal,  S  is  a  diagonal  real  k  x  k  matrix  whose  entries  are  all 
nonnegative,  and  Ck+i  is  the  {k  +  1)*^*  greatest  singular  value  of  A.  To  do  so,  we  first  use 
the  algorithm  of  Subsection  4.1  to  construct  the  matrices  B  and  P  in  (84),  (91),  and  (93). 
Then,  we  make  the  following  four  steps: 

1.  Construct  a  lower  triangular  real  k  x  k  matrix  L,  and  a  real  n  x  k  matrix  Q  whose 
columns  are  orthonormal,  such  that 

P  =  LQ"^  (98) 

(see  Remark  18  for  details  concerning  the  construction  of  such  matrices  L  and  Q). 
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Form  the  product 

C  =  BL. 

(99) 

Construct  an  SVD  of  C,  that  is. 

C  =  U^W^, 

(100) 

where  U  is  a  real  m  x  k  matrix  whose  columns  are  orthonormal,  S  is  a  diagonal 
k  X  k  matrix  whose  entries  are  all  nonnegative,  and  W  is  a  real  k  x  k  matrix  whose 
columns  are  orthonormal  (see,  for  example.  Chapter  8  in  [15]  for  details  concerning  the 
construction  of  such  an  SVD). 

4.  Form  the  product 

V  =  QW.  (101) 

The  matrices  U,  S,  and  V  obtained  via  the  preceding  four  steps  (after  hrst  using  the  al¬ 
gorithm  of  Subsection  4.1  to  construct  the  matrices  B  and  P)  satisfy  (97);  see  Remark  23 
below. 

Remark  22  Steps  2  and  4  in  the  procednre  of  the  present  snbsection  are  somewhat  snbtle 
nnmerically.  Both  Steps  2  and  4  involve  constructing  products  of  matrices,  and  in  general 
constructing  the  product  S'  T  of  matrices  S  and  T  can  be  numerically  unstable.  Indeed,  in 
general  some  entries  of  S'  or  T  can  have  unmanageably  large  absolute  values,  while  in  exact 
arithmetic  no  entry  of  the  prodnct  S'  T  has  an  unmanageably  large  absolute  value;  in  snch 
circnmstances,  constructing  the  prodnct  S'  T  can  be  unstable  in  hnite-precision  arithmetic. 
However,  this  problem  does  not  arise  in  Steps  2  and  4  above,  dne  to  (84),  (81),  the  fact  that 
the  columns  of  B  constitnte  a  snbset  of  the  colnmns  of  A  (so  that  ||i?||  <  ||R||),  and  the  fact 
that  the  columns  of  Q  are  orthonormal,  as  are  the  columns  of  W. 

Remark  23  In  this  remark,  we  demonstrate  that  the  matrices  U,  S,  and  V  satisfy  (97). 
Indeed,  suppose  that  f3  and  7  are  positive  real  nnmbers  such  that  7  >  1  and  y  dehned  in  (90) 
is  nonnegative.  Then,  combining  (91),  (75),  (98),  (99),  (100),  (101),  and  (79)  yields 

\\UT.V^  -A\\ 

<  7^  -I-  1  ^a/4/c  {n  —  k)  +  1  -|-  1  j  -1-/37  \/2lm  ak+i  (102) 

with  probability  not  less  than  y  dehned  in  (90),  where  ak+i  is  the  {k  +  1)®*  greatest  singular 
value  of  A.  The  bound  (102)  is  a  precise  version  of  (97).  As  in  Subsection  4.1,  if  we  choose 
/3  =  3/4,  7^  =  5,  and  I  =  k  +  20,  and  combine  (102)  and  (90),  then  we  obtain 

\\UYi  —  A||  <  10  A/Ar(Ar+~20yr?m  ak+i  (103) 

with  probability  not  less  than  1  —  10“^^.  Table  1  contains  similar  results  obtained  by  taking 
other  valnes  for  I  —  k,  {3,  and  7. 
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Remark  24  When  the  singular  values  of  A  decay  sufficiently  fast,  and  I  is  much  less  than 
m,  the  factors  ^/2lm(3‘^  7^  +  1  and  y/2lm  in  (102)  are  much  larger  than  necessary.  Indeed, 
suppose  that  j  is  a  positive  integer  with  k  +  j  <  m  and  k  +  j  <  n,  and  f3  and  7  are  positive 
real  numbers,  such  that  7  >  1  and  <h  +  4'  >  1,  where  <h  is  dehned  in  (39),  and  is  dehned 
in  (64).  Then,  combining  (93),  (75),  (98),  (99),  (100),  (101),  and  (79)  yields 

\\UT.V^  -  A\\  <iak+i+riak+j+i  (104) 

with  probability  not  less  than  $  +  T  —  1,  where  <h  is  dehned  in  (39),  T  is  dehned  in  (64), 
(Tfc+i  is  the  {k  +  1)®*  greatest  singular  value  of  A,  au+j+i  is  the  {k  +  j  +  1)®*  greatest  singular 
value  of  A,  ^  is  dehned  in  (94),  and  rj  is  dehned  in  (95).  When  j,  k,  and  I  are  all  much  less 
than  m,  clearly  ^  is  much  less  than  7. 

4.3  CPU  time  and  memory  requirements 

In  this  subsection,  we  tabulate  the  numbers  of  hoating-point  operations  and  words  of  memory 
required  by  the  algorithms  described  in  Subsections  4.1  and  4.2,  as  applied  once  to  a  matrix 

4. 

4.3.1  Interpolative  decomposition 

The  algorithm  described  in  Subsection  4.1  incurs  the  following  costs  in  order  to  compute  an 
approximation  to  an  interpolative  decomposition  of  A: 

1.  Forming  R  in  (86)  requires  applying  to  I  vectors. 

2.  Computing  S  and  P  in  (87)  or  (96)  costs  0{lkn  log(/)). 

3.  Forming  B  in  (88)  requires  applying  A  to  k  vectors,  where  each  vector  has  a  single 
entry  of  1  and  n  —  1  entries  of  0. 

Summing  up  the  costs  in  Steps  1-3  above,  we  conclude  that  the  algorithm  of  Subsection  4.1 
costs 

Cm  =  A;  ■  Ca  +  /  ■  +  0{lkn  log(/)),  (105) 

where  Ca  is  the  cost  of  applying  4  to  a  real  n  x  1  column  vector,  and  C^t  is  the  cost  of 
applying  A"^  to  a  real  m  x  1  column  vector. 

4.3.2  Singular  Value  Decomposition 

The  algorithm  described  in  Subsection  4.2  incurs  the  following  costs  in  order  to  compute 
an  approximation  to  an  SVD  of  A  from  the  matrices  B  and  P  obtained  via  the  algorithm 
described  in  Subsection  4.1: 

1.  Computing  the  matrices  L  and  Q  in  (98)  costs  0{k‘^n). 

2.  Forming  C  in  (99)  costs  0{k‘^m). 

3.  Computing  the  SVD  (100)  of  C  costs  0{k‘^m). 
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4.  Forming  V  defined  in  (101)  costs  0{k‘^n). 

Summing  up  the  costs  in  Steps  1-4  above,  we  conclude  that  the  algorithm  of  Subsection  4.2 
costs 

CsvD  =  C'lD  +  {m  +  n)),  (106) 

where  Cm  is  the  cost  in  (105)  of  applying  the  algorithm  described  in  Subsection  4.1  in  order 
to  obtain  B  and  P. 

Remark  25  We  observe  that  in  order  to  form  the  matrix  G  Am  (86),  we  need  only  to  take 
each  column  of  A  and  form  its  inner  products  with  the  I  rows  of  G.  In  order  to  form  B 
in  (88),  we  need  only  to  collect  together  certain  columns  of  A.  Thus,  we  only  ever  need  to 
access  columns  of  A,  and  never  need  to  transpose  A  explicitly.  The  data  movement  required 
by  the  algorithm  of  the  present  paper  is  hence  fairly  simple. 

4.4  Description  of  Table  1 

Tables  1. 1-1.6  provide  an  upper  bound  on  the  probability  that 

\\UEV^  -  A\\  >  (Vkh^ak+i,  (107) 

where  U,  S,  and  V  are  the  matrices  in  the  approximation  to  an  SVD  of  the  mx  n  matrix  A 
in  (102).  In  (107),  k  and  I  are  any  positive  integers  with  k  <  I,  such  that  I  <  m  and  I  <  n, 
(Tfc+i  is  the  {k  +  l)*^*  greatest  singular  value  of  A,  and  (  takes  on  the  values  specihed  by  the 
penultimate  columns  of  the  tables.  The  quantity  is  dehned  by  the  formula 

^l-k,P,y  =  1  —  X)  (108) 

where  y  is  dehned  in  (90),  and  I  —  k,  P,  and  7  take  on  the  values  specihed  by  the  hrst, 
second,  and  third  columns  of  the  tables.  The  quantity  (  is  specihed  by  P  and  7  via  (102). 
Please  note  that  If^.^,^^  depends  only  on  /  —  /c,  P,  and  7,  and  provides  an  upper  bound 
that  is  otherwise  independent  of  k,  m,  n,  and  A]  (103)  provides  a  similar  bound.  When  the 
singular  values  of  A  decay  sufficiently  fast,  and  I  is  much  less  than  m,  the  factor  of  ^/m  in 
the  right-hand  side  of  (107)  is  larger  than  necessary;  see  Remark  24  above. 

Remark  26  Due  to  (91),  the  quantity  dehned  in  (108)  also  provides  an  upper  bound 

on  the  probability  that 

\\B  P  —  A\\  >  (  Vklnm  ak+i,  (109) 

where  B  and  P  are  the  matrices  in  the  approximation  to  an  interpolative  decomposition  of 
the  mx  n  matrix  A  in  (91). 

5  Numerical  results 

In  this  section,  we  describe  the  results  of  four  numerical  tests  of  the  algorithm  of  the  present 
paper.  Table  2  summarizes  the  numerical  output  of  the  examples  described  in  the  present 
section. 
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Tables  2. 1-2.4  display  the  results  of  applying  the  algorithm  of  the  present  paper  once 
to  a  real  n  x  n  matrix  A,  for  the  indicated  values  of  n.  The  matrix  A  is  dehned  at  the 
end  of  the  present  section.  The  numbers  k  and  /  are  those  from  Section  4;  k  is  the  rank 
of  the  approximations  to  A,  and  /  is  the  number  of  rows  in  the  matrix  G  whose  entries 
are  i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance  (the  algorithm  uses  the 
product  G  A).  The  displayed  times  refer  to  the  seconds  of  CPU  time  used  by  the  algorithm  to 
compute  both  the  approximation  to  an  interpolative  decomposition  and  the  approximation 
to  an  SVD  of  A.  (Please  note  that  our  implementation  is  optimized  for  accuracy  and  for 
analyzing  the  numerical  properties  of  the  algorithm,  and  is  probably  not  very  efficient.)  The 
numbers  and  Uk+i  are  those  from  the  dehnition  of  A  below;  furthermore,  appears  in 
the  bounds  (91)  and  (102)  on  the  errors  of  the  approximations.  The  number  (5id  is  the  norm 
of  the  difference  between  A  and  an  approximation  5  P  to  an  interpolative  decomposition  of 
A,  that  is, 

=  \\B  P  -  Al  (110) 

where  the  matrices  B  and  P  are  those  from  (91).  The  number  (5svd  is  the  norm  of  the 
difference  between  A  and  the  approximation  U  S  to  an  SVD  of  A,  that  is, 

5svd  =  lies  All,  (111) 

where  the  matrices  U,  S,  and  V  are  those  from  (102). 

We  dehne  (5rei.  max.  as  follows.  First,  we  dehne  auxiliary  vectors  A,  . . . ,  P,  and 
r^,  r^,  . . . ,  t\  where  j  =  30  in  Examples  1  and  2,  j  =  110  in  Example  3,  and  j  =  6  in 
Example  4.  We  choose  the  test  vectors  . . . ,  P~^,  P  to  include  a  variety  of  deterministic 

and  random  vectors  (specihcally,  we  set  every  entry  of  A  to  be  1,  and  use  a  random  number 
generator  to  generate  ...,  P~^,  P  so  that  their  entries  are  i.i.d.  Gaussian  random 

variables  of  zero  mean  and  unit  variance).  For  any  i  =  1,  2,  . . . ,  j  —  1,  j ,  we  dehne  r*  to  be 
the  vector  resulting  from  the  application  of  P  P  —  A  to  P,  that  is, 

P  =  {BP-A)t\  (112) 


where  the  matrices  B  and  P  are  those  from  (91).  Then,  we  dehne  (5rei.max.  via  the  formula 


<^rel. 


max. 


max 


/ maXp=i,2,...,m-i,m  I('r0p|\ 

V  maXg=i,2,...,n-l,n  \{P)q\  J  ’ 


(113) 


where  {P)p  is  the  entry  of  r*,  and  {P)q  is  the  entry  of  P. 

All  estimates  displayed  in  Table  2  are  the  maximum  values  obtained  from  three  indepen¬ 
dent  realizations  of  the  random  variables  involved. 

The  values  of  5id  and  (5svd  displayed  in  Tables  2.2,  2.3,  and  2.4  are  those  obtained  via  the 
power  method  for  estimating  the  norm  of  a  matrix,  after  the  estimates  stabilized  to  three 
signiheant  hgures.  The  values  of  (5id  and  (5svd  displayed  in  Table  2.1  are  those  obtained  after 
100  iterations  of  the  power  method.  The  estimates  of  (5id  and  (5svd  summarized  in  Table  2.1 
did  not  stabilize  to  three  signiheant  hgures  after  100  (or  any  other  number  of)  iterations, 
undoubtedly  due  to  round-oh. 

We  performed  all  computations  using  IEEE  standard  double-precision  variables,  whose 
mantissas  have  approximately  one  bit  of  precision  more  than  16  digits  (so  that  the  relative 
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precision  of  the  variables  is  approximately  .2e-15).  We  ran  all  computations  on  a  2.8  GHz 
Pentium  Xeon  microprocessor  with  512  KB  of  L2  cache  and  2  GB  of  RAM.  We  compiled  the 
Fortran  77  code  using  the  Lahey-Fujitsu  compiler,  with  the  optimization  flag  — o2  enabled. 

In  our  implementation,  we  computed  SVDs  using  2-sided  plane  (Jacobi/Givens)  rotations 
(see,  for  example,  Ghapter  8  in  [15]).  We  used  an  algorithm  based  upon  pivoted  “QR”  de¬ 
compositions  to  compute  the  matrices  S  and  P  in  (87)  and  (96)  (see,  for  example,  Ghapter  5 
in  [15]  for  a  description  of  R”  decompositions,  and  [6]  for  further  details  regarding  our 
particular  implementation). 

In  Examples  1,  2,  and  3,  we  use  a  pseudorandom  number  generator  to  construct  real 
n  X  1  vectors  jR,  . . . ,  /ib  and  z/^,  z/^,  . . . ,  z/-^"^,  z/-^,  such  that  their  entries  are  a 

realization  of  i.i.d.  Gaussian  random  variables  of  zero  mean  and  unit  variance,  with  j  =  20 
in  Examples  1  and  2,  and  j  =  60  in  Example  3.  We  orthonormalize  /i^,  /i^,  . . . , 
via  the  Gram-Schmidt  process  with  reorthogonalization  (see,  for  example,  [3])  to  obtain  real 
n  X  1  vectors  . . . ,  u\  and  do  the  same  with  z/^,  z/^,  . . . ,  z/-^“^,  z/-^  to  obtain 

real  n  x  1  vectors  v^,  . . . ,  v^~^,  vK  We  denote  by  ai,  <72,  . . . ,  (Jj-i,  Uj  the  positive  real 

numbers  displayed  in  Figure  1  (we  use  the  numbers  in  Figure  1.1  for  Example  1,  the  numbers 
in  Figure  1.2  for  Example  2,  and  those  in  Figure  1.3  for  Example  3).  We  dehne  A  to  be  the 
n  X  n  matrix  given  by  the  formula 


j 

A  =  ^u^  a,  {Pf. 

i=l 


(114) 


Glearly,  the  rank  of  A  is  j.  Since  P,  v?,  . . . ,  P  are  orthonormal,  as  are  v‘^,  . . . , 

v^~^,  ,  the  singular  value  of  A  is  (Tj,  for  alH  =  1,  2,  . . . ,  j  —  1,  j.  Table  2  displays  the 

results  of  applying  the  algorithm  of  the  present  paper  to  A,  for  various  values  of  n  (Table  2.1 
displays  the  results  for  Example  1,  Table  2.2  displays  the  results  for  Example  2,  and  Table  2.3 
displays  those  for  Example  3). 

Example  4  is  designed  to  illustrate  that  factors  of  the  order  of  ^J4:k{n  —  k)  +  1  are 
necessary  in  bounds  such  as  (91),  (93),  (102),  and  (104).  This  example  is  identical  to 
Examples  1,  2,  and  3,  using  the  same  matrix  A  defined  in  (114),  but  with  n  assumed  to  be 
divisible  by  8,  with  j  =  4,  and  using  the  numbers  ui,  (T2,  <73,  and  <74  displayed  in  Figure  1.4 
((7i  =  1,  (72  =  1,  (73  =  .le-7,  and  (74  =  .le-7),  and  with  the  following  vectors  vA,  u^, 
and  v^,  v^, 

1  1  ...  1  1  ),  (115) 

y  n 

=  -1  1  -1  ...  1  -1  1  -1),  (116) 
\/n 


(„T  =  ^(i  1-1-111-1  -1  ...  11-1-111-1  -1). 

(117) 

(„T  =  ^(i  1  1  1  -1  -1  -1  -1  ...  1  1  1  1  -1  -1  -1  -1). 

(118) 

=  1  1  ...  1  1  1  0  ),  (119) 
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(i;2)T  =  (0  0  ...  0  0  0  1), 

_1  1  _1  ...  1-11-10  0), 

{v^f  =  ^  (  1  0  -1  0  0  0  . . .  0  0  )  , 

that  is, 


(120) 

(121) 

(122) 


1.  (a)  every  entry  of  is  1/^/n, 

(b)  entries  1,  2,  . . . ,  n  —  2,  n  —  1  of  are  1/ ^/n  —  1,  and  entry  n  of  is  0, 

2.  (a)  every  even  entry  of  is  — l/y^,  and  every  odd  entry  of  v?  is  1/^/n, 

(b)  entries  1,  2,  . . . ,  n  —  2,  n  —  1  of  are  0,  and  entry  n  of  v"^  is  1, 

3.  (a)  the  hrst  pair  of  entries  of  v?  is  1/^/n,  the  second  pair  of  entries  is  —1/^/n,  the 

third  pair  of  entries  is  l/v^i  fhe  fourth  pair  of  entries  is  —  1/y^,  and  so  on  (with 
each  successive  pair  of  entries  alternating  sign), 

(b)  every  even  entry  of  except  for  entry  n  is  —l/\/n  —  2,  every  odd  entry  of 
except  for  entry  n  —  1  is  l/\/n  —  2,  and  entries  n  —  1  and  n  of  are  0, 

4.  (a)  the  hrst  quadruplet  of  entries  of  is  l/y/n,  the  second  quadruplet  of  entries  is 

— l/yhr,  the  third  quadruplet  of  entries  is  1/v^i  fhe  fourth  quadruplet  of  entries 
is  —  l/ydr,  and  so  on  (with  each  successive  quadruplet  of  entries  alternating  sign), 
and 

(b)  entry  1  of  is  1/v^,  entry  2  of  is  0,  entry  3  is  — 1/v^,  and  entries  4,  5,  ... , 
n  —  1,  n  are  0. 


For  this  example  (Example  4), 


(Ti  (n^)’’"  +  (72  (n^) 


2\T 


1 

^/n  {n  -  1) 


/I  1 
1  1 
1  1 
1  1 


\  1  1 


1  1  \ 
1  1  -Vn^ 

1  1  1 
1  1  -^/n^ 

1  1 

1  1  -y/n^  / 


(123) 


Clearly,  v?,  u^,  and  u‘^  are  orthonormal,  as  are  v^,  v^,  and  Therefore,  the 

singular  value  of  A  is  o-j,  for  all  i  =  1,  2,  3,  4.  Table  2.4  displays  the  results  of  applying  the 
algorithm  of  the  present  paper  to  A,  for  various  values  of  n. 


Remark  27  All  numerical  data  that  we  have  examined  —  including  the  data  displayed  in 
Table  2,  as  well  as  the  data  from  further  experiments  —  appear  to  satisfy  the  bounds  (91), 
(93),  (102),  and  (104).  However,  we  are  not  entirely  certain  that  bounds  such  as  (91),  (93), 
(102),  and  (104)  must  necessarily  contain  factors  of  the  order  of  ^/m  and 
We  would  not  be  surprised  if  tighter  bounds  were  possible. 
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The  loss  of  precision  displayed  in  Table  2.1  as  n  increases  is  probably  largely  due  to  round¬ 
off  (compare  Table  2.2),  whereas  the  loss  of  precision  displayed  in  Table  2.4  as  n  increases 
suggests  that  any  bounds  such  as  (91),  (93),  (102),  and  (104)  must  contain  a  factor  of  the 
order  of  In  contrast,  in  many  practical  situations  the  bounds  mentioned 

in  Remarks  20  and  24  are  effectively  independent  of  m. 
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10  3  10^ 

88  4  W 

790  5  W 

6,600  7  KF 

58,000  9  W 


.15e-0 

.18e-l 

.21e-2 

.18e-3 

.19e-4 


p  I  7  II 

7.9  5  I  10^  .85e-2 

72  .lle-3 

620  .15e-5 

5,500  KTFF  .17e-7 

53,000  11  FO®  .19e-9 


2 

7.9 

5 

o 

to 

2 

66 

7 

10^ 

2 

589 

9 

10^ 

2 

5,300 

11 

10^ 

2 

51,000 

12 

10*^ 

1.02 

3 

10 

.18e-0 

0.88 

4 

10 

.16e-l 

0.88 

4 

10 

.lOe-4 

0.79 

5 

10 

.16e-12 

0.72 

6 

10 

.59e-32 

0 

10 

3 

100 

1 

7.9 

5 

100 

2 

7.9 

5 

100 

4 

7.2 

6 

100 

8 

6.6 

7 

100 

16 

6.6 

7 

100 

32 

6.2 

8 

100 

10 

3 

1,000 

.18e-l 

72 

6 

1,000 

.lle-3 

66 

7 

1,000 

.61e-6 

62 

8 

1,000 

.93e-ll 

58 

9 

1,000 

.37e-21 

55 

10 

1,000 

.73e-44 

53 

11 

1,000 

.14e-93 

Table  1  (See  Subsection  4.4.) 


k 

1 

n 

time  (sec.) 

O'k 

o'fc+i 

^rel.  max. 

<5id 

<^SVD 

10 

10 

10^ 

.29e-03 

.232e-15 

.200e-15 

.162e-14 

.533e-14 

.534e-14 

10 

10 

10^ 

.27e-02 

.232e-15 

.200e-15 

.211e-14 

.286e-14 

.934e-14 

10 

10 

10^ 

.42e-01 

.232e-15 

.200e-15 

.844e-14 

.187e-13 

.187e-13 

10 

10 

10® 

.16e+01 

.232e-15 

.200e-15 

.181e-13 

.530e-13 

.189e-12 

10 

10 

10® 

.12e+02 

.232e-15 

.200e-15 

.856e-13 

.121e-12 

.683e-ll 

(2.1) 


k 

1 

n 

time  (sec.) 

O'k 

O'k+l 

^rel.  max. 

(^ID 

<^SVD 

10 

10 

10^ 

.28e-03 

.152e-07 

.lOOe-07 

.139e-06 

.156e-06 

.156e-06 

10 

10 

10® 

.28e-02 

.152e-07 

.lOOe-07 

.953e-07 

.806e-07 

.806e-07 

10 

10 

10^ 

.40e-01 

.152e-07 

.lOOe-07 

.104e-06 

.242e-06 

.242e-06 

10 

10 

10® 

.16e+01 

.152e-07 

.lOOe-07 

.165e-06 

.182e-06 

.182e-06 

10 

10 

10® 

.12e+02 

.152e-07 

.lOOe-07 

.102e-06 

.134e-06 

.134e-06 

(2.2) 


k 

/ 

n 

time  (sec.) 

<^k 

o'fc+i 

^rel.  max. 

<^ID 

<^SVD 

30 

30 

50,000 

.41e+01 

.206e-08 

.lOOe-08 

.179e-07 

.482e-07 

.482e-07 

30 

31 

50,000 

.42e+01 

.206e-08 

.lOOe-08 

.213e-07 

.348e-07 

.348e-07 

30 

32 

50,000 

.41e+01 

.206e-08 

.lOOe-08 

.114e-07 

.304e-07 

.304e-07 

30 

34 

50,000 

.43e+01 

.206e-08 

.lOOe-08 

.139e-07 

.201e-07 

.201e-07 

30 

38 

50,000 

.45e+01 

.206e-08 

.lOOe-08 

.104e-07 

.144e-07 

.144e-07 

30 

46 

50,000 

.48e+01 

.206e-08 

.lOOe-08 

.805e-08 

.863e-08 

.868e-08 

k 

1 

n 

time  (sec.) 

(2. 

3) 

O'k+l 

^rel.  max. 

f^ID 

5s  VD 

2 

2 

.48  ■  10^ 

.19e-04 

.lOOe+Ol 

.lOOe-07 

.128e-06 

.757e-07 

.757e-07 

2 

2 

oo 

o 

CO 

.13e-03 

.lOOe+01 

.lOOe-07 

.374e-06 

.260e-06 

.260e-06 

2 

2 

OO 

o 

.22e-02 

.lOOe+01 

.lOOe-07 

.728e-05 

.663e-05 

.663e-05 

2 

2 

.48  ■  10® 

.16e-01 

.lOOe+01 

.lOOe-07 

.853e-05 

.714e-05 

.714e-05 

2 

2 

.48  ■  10® 

.18e-00 

.lOOe+01 

.lOOe-07 

.152e-04 

.114e-04 

.114e-04 

2 

2 

OO 

o 

-1 

.28e+01 

.lOOe+01 

.lOOe-07 

.306e-04 

.216e-04 

.216e-04 

(2.4) 


Table  2  (See  Section  5.) 
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